
/*************
author:: libo
license:: MIT
*****************/
#include<Math/LB_Math.h>
#define factorial LB_factorial


// void matrix_multiply(double * __restrict A, double * __restrict B, double * __restrict C, int N) {
//     for (int i = 0; i < N; i++) {
//         for (int j = 0; j < N; j++) {
//             double sum = 0.0;
//             // 注意：这里没有使用 OpenMP 指令
//             for (int k = 0; k < N; k++) {
//                 sum += A[i * N + k] * B[k * N + j];
//             }
//             C[i * N + j] = sum;
//         }
//     }
// }

// // 内存分配时确保数组以适当边界对齐
// double* allocate_aligned_matrix(int N, size_t alignment) {
// #ifdef _WIN32
//     return (double*)_aligned_malloc(N * N * sizeof(double), alignment);
// #else
//     return (double*)aligned_alloc(alignment, N * N * sizeof(double)); // 对于非Windows平台
// #endif
// }

// // 释放对齐的内存
// void free_aligned_matrix(double* matrix) {
// #ifdef _WIN32
//     _aligned_free(matrix);
// #else
//     free(matrix); // 对于非Windows平台，如果使用 aligned_alloc，则可以直接用 free
// #endif
// }
int factorial(int n)
{
    int re=1;
    for(int i=2;i<=n;i++)
    {
        re*=i;
    
    }
    return re;

}
int *pailie_map_i_to_reverse_order(int n)
{
    int num=factorial(n);
    int *re=(int*)malloc(sizeof(int)*n*num);

    for(int i=0;i<num;i++)
    {
        int chu=i;
        for(int j=0;j<n;j++)
        {
            int yu=chu%(j+1);
            chu=chu/(j+1);
            re[i*n+n-1-j]=yu;
        }
    }
    return re;
}
//Convert the result of the above function


int * map_reverse_order_to_pailie(int* r,int n)
{
    int num=factorial(n);
    int *re=(int*)malloc(sizeof(int)*n*num);
    for(int i=0;i<num;i++)
    {
        for(int j=0;j<n;j++)
        { 
            re[i*n+j]=j;
        } 
    }
    int ix,iy,temp;
    for(int i=0;i<num;i++)
    {
        for(int j=0;j<n;j++)
        {
            ix=i*n+j;
            iy=i*n+j+r[i*n+j];
            
            for(int l=iy;l>ix;l--)
            { 
                temp=re[l-1];
                re[l-1]=re[l];
                re[l]=temp;
            }          
        }   
    }
    return re;
}
//compute the station pick m elements from n elements

static int increase_a_step_for_combination(int m,int n,int* temp_array)
{
	int temp_i;
    int *mu_array=temp_array;
	for(int i=0;i<m;i++)
	{
		temp_i=mu_array[i]+1;
        	if(temp_i<n)
        	{
            		if((i+1)<m)
            		{
                		if(temp_i>=mu_array[i+1])
                		{
                    			continue;
                		}
            		}
            
            		int* temp_mu=&(mu_array[m]);
            		memmove(temp_mu, mu_array,sizeof(int)*m);
            		temp_mu[i]=temp_i;
            		mu_array=temp_mu;
            		for(int j=0;j<i;j++)
            		{
                		mu_array[j]=j;
            		}
/*Execute once and return*/
            		return 1;

        	}
        	else
        	{
            		return 0;
        	}
	}
    return 0;
    
}
int* compute_combination(int m,int n,int* re,int* re_len)
{
    if(m>n){return re;}
    if(re==NULL)
    {
        *re_len=1;
        for(int i=0;i<m;i++)
        {
            (*re_len)*=(n-i);     
        }
        (*re_len)=(*re_len)/factorial(m);
    //printf("%d\n",len);
        re=(int*)malloc(sizeof(int)*(*re_len)*m);
    }
    for(int i=0;i<m;i++)
    {
        re[i]=i;
    } 
    for(int i=0;i<(*re_len)-1;i++)
    {
        increase_a_step_for_combination(m,n,&(re[i*m]));
    }

    return re;
}
// int * compute_combination(int m,int n)
// { 
// 	int * re=0;
// 	if(m>n)
// 	{return re;}
// 	int len=1;
// 	for(int i=0;i<m;i++)
// 	{
// 		len*=(n-i);		

// 	}
// 	len=len/factorial(m);
// 	//printf("%d\n",len);
// 	re=(int*)malloc(sizeof(int)*len*m);
    
//     for(int i=0;i<m;i++)
//     {
//         re[i]=i;
//     } 
//     for(int i=0;i<len-1;i++)
//     {
    
//         increase_a_step_for_combination(m,n,&(re[i*m]));
//     }
// 	return re;
// }


double inner_product(double* a,double* b,int len )
{   
    double re=0;
    for(int i=0;i<len;i++)
    {
        re+=a[i]*b[i];
    }
    return re;
}

int normalize(double *a,int len)
{
    double norm =0;
    for(int i=0;i<len;i++)
    {
        norm+=a[i]*a[i];
    }
    norm=SAFE_SQRT(norm);
    //printf("norma:%.9lf\n",norm );
    if(double_equal_zero(norm))
    {
        return 0;
    }
    for(int i=0;i<len;i++)
    {
        a[i]/=norm;
    }
    return 1;
}
// 
// min 0.5xTHx+cTx
//  Ax<=b
//  Ex=d
// 
void compute_convex_optimize(double**H,double*c,
    double**A,double* b,double **E,double *d,
    double*x,int len)
{

}



#undef factorial
